SPECTRUM OF THE JACOBI TAU APPROXIMATION FOR THE 
SECOND DERIVATIVE OPERATOR 
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Abstract. It is proved that the eigenvalues of the Jacobi Tau method for the second derivative 
operator with Dirichlet boundary conditions are real, negative and distinct for a range of the Jacobi 
parameters. Special emphasis is placed on the symmetric case of the Gegenbauer Tau method where 
the range of parameters included in the theorems can be extended and characteristic polynomials 
given by successive order approximations interlace. This includes the common Chebyshev and Leg- 
endre, Tau and Galerkin methods. The characteristic polynomials for the Gegenbauer Tau method 
are shown to obey three term recurrences plus a constant term which vanishes for the Legendre Tau 
and Galerkin cases. These recurrences are equivalent to a tridiagonal plus one row matrix structure. 
The spectral integration formulation of the Gegenbauer Tau method is shown to lead directly to 
that fundamental and well-conditioned tridiagonal plus one row matrix structure. A Matlab code is 
provided. 
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1. Introduction. Constructing polynomial approximations to solutions of dif- 
ferential equations is the basic ingredient of most numerical methods. Approximations 
based on orthogonal polynomials have been widely used (e.g. [3], because their 
rate of convergence is faster than algebraic for arbitrary boundary conditions when 
the solution is smooth. The purpose of this paper is to give a rigorous proof that 
the spectrum of the Jacobi Tau approximation is real, negative and distinct for the 
second order operator with Dirichlet boundary conditions. The Jacobi Tau class of 
spectral methods includes the common Chebyshev and Legendre Tau and Galerkin 
formulations, as demonstrated below. The general method of proof is similar to that 
used by Gottlieb and Lustman |5J[7] to prove such results for the Chebyshev colloca- 
tion operator. However, we argue in section ETT1 that Gottlieb and Lustman's proof 
for the collocation operator is not complete. 

The spectrum of Jacobi Tau approximation for the 1st order operator has been 
considered elsewhere 0]. Here, we consider polynomial approximations to the eigen- 
value problem 

(1.1) < L^-^\u -Kx<l with u(±l) = 0. 

dx 

The spectrum of Jacobi polynomial approximations to this eigenvalue problem is 
directly relevant to numerical simulations of the diffusion equation u t = u xx which is 
itself a building block for numerical solution of various other problems including the 
Stokes and Navier-Stokes equations (e.g. J3J §5.1], |18|). 

The 2nd order problem 1)1. ljl is a self-adjoint, negative definite Sturm-Liouville 
differential eigenproblem, so its eigenvalues A are real, negative and distinct. The 
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eigenmodes separate into even and odd modes and have the simple exact expressions 
u e (x) = cos(2fc — A 

(1.2) 

u (x) = sin kirx, A 

for k= 1,2,3,... . 

If u n {x) is a polynomial approximation of degree n to the exact solution u(x), 
then u n {x) satisfies the following differential equation 

(1.3) \u n {x) - D 2 u n {x) = Unix) 

where the residual R n (x) is a polynomial of degree n in x and D = d/dx. We can 
invert this relation to express the polynomial approximation u n (x) in terms of the 
residual R n (x) 

[n/2] 

(1.4) u n {x)= t iY jt i k D 2k R n {x) 

k=Q 

where ^ = 1/A, and [n/2] denotes the greatest integer less or equal to n/2. We 
can assume that A ^ because A = with u„(±l) = necessarily corresponds to 
the trivial solution U n [x) = 0, Vx in [—1,1], as shown below. The inversion (|1.4fl 
follows from formal application of the geometric (Neumann) series for (1 — /iD 2 ) -1 = 
Sfco ^ k D 2k which terminates since Rn{x) is a polynomial. That inversion can also 
be derived by repeated application of the operator [iD 2 to equation 1)1. 3JI . Summation 
of the resulting suite of equations leads to (|1.4)l thanks to telescopic cancelations on 
the left hand side. 

Spectral methods fit in the general framework of the method of weighted residuals 
In the Tau method §10.4.2], the polynomial approximation u n (x) is deter- 
mined from the boundary conditions u„(±l) = and the requirement that R n (x) is 
orthogonal to all polynomials p n ~2{x) of degree n — 2 or less with respect to a weight 
function W(x) > in the interval (—1, 1) 

(1.5) J R n (x) Pn - 2 (x)W(x)dx = 0. 

These requirements provide n + 1 equations for the n + 1 undetermined constants 
in the polynomial approximation u n (x). For the Jacobi weight function W a ^{x) = 
(1 - x) a {l + xf, the residual (Q can be written as 

(1.6) R n (x) = T XP^\x) + Tl XPi a Jp(x) 

for some ^-independent coefficients To and n, where Pn a '^\x) is the Jacobi polynomial 
of degree n fsect. |BT|) . This follows from orthogonality of the Jacobi polynomials in 
— 1 < x < 1 with respect to the Jacobi weight W ay p(x) which implies orthogonality of 
the Jacobi polynomial of degree k to any polynomial of degree k— 1 or less with respect 
to that weight function. Jacobi polynomials are the most general class of polynomial 
solutions of a Sturm-Liouville eigenproblcm that is singular at ±1 as required for 
faster than algebraic convergence §9.2.2, §9.6.1]. It is now easy to verify from 
(tOl) and (JUl) that if A = then D 2 u n (x) = R n (x) = for all x in (-1,1) but 



= -(2k-l) 2 '-, 

= -k 2 7T 2 . 
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the boundary conditions u„(±l) = would then require that u n (x) = for all x in 
[—1, 1]. Therefore we can assume that A 7^ 0. 

In the Galerkin approach, u n {x) is determined from the boundary conditions 
u n (±l) = and orthogonality of the residual R n (x) to all polynomials of degree n 
that vanish at x = ±1, with respect to a weight function W(x) > 0. In other words, 
the test functions are in the same space (polynomials of degree n) as the trial functions 
and they satisfy the same boundary conditions. Such polynomials can be written in 
the form (1 — x 2 )p n -2(x) where p n -2(x) is an arbitrary polynomial of degree n — 2, 
and the Galerkin equations can be written as 

(1.7) J R n (x){l-x 2 )p,^ 2 (x)W{x)dx = 0. 

For the Jacobi weight W(x) = W a ,^{x) = (1 — x) a (l + x) 13 , the Galerkin method is 
therefore equivalent to the Tau method for the weight W a +i t f)+i_(x) and the residual 
controlled by (|1 . T|) has the form 

(1.8) Rn(x) = T XP^ +1 \x)+T 1 XPi a _ + 1 ^ +1 \x). 

This residual can be written in terms of the derivatives of pj^ 1 '^ +1 \x) and Pn''^\x) 
by making use of i|B.7f) . Since we consider a range of parameters a and (3, the Jacobi- 
Tau method also includes some Jacobi-Galerkin methods. 

In the collocation approach, u n (x) is determined from the boundary conditions 
it n (±l) = and enforcing R n (xj) = at the n — 1 interior Gauss-Lobatto points Xj 
such that DP^ a,P) {xj) = 0, j = 1, . . . , n - 1 3, §2.2]. The residual JOJ takes the 
form eqn. (4.5)] 

(1.9) R n (x) = (A + Bx)DP^\x), 

for some A and B independent of x. The collocation residual i|1.9|) is provided for 
completeness since we do not have results about the collocation method and we raise 
doubts about the validity of the proof proposed in [?]■ That residual can be written 
in several equivalent forms by using the properties of Jacobi polynomials (sect. iBr]) . 

The characteristic polynomials for the eigenvalues \x — 1/A are derived in section 
121 from the explicit expression (|1.4f> for u n (x) in terms of the residual R n (x) whose 
form is specified by the Jacobi Tau (or Galerkin) method as in (|1.6I) and l|1.8|) . The 
zeros of these characteristic polynomials are shown to be real, negative and distinct 
in section |3J Recurrence relations for the Gegenbauer Tau characteristic polynomials 
are derived in section 0] where it is shown that the underlying fundamental matrix 
structure is tridiagonal + one row. Section[S]discusses some implementation issues and 
shows that the spectral integration implementation directly leads to the tridiagonal + 
one row structure which is well-conditioned. Some of the key properties of Jacobi and 
Gegenbauer polynomials used in this paper are summarized in appendix^] We use a 
non-standard normalization for Gegenbauer polynomials, denoted Gn\x), since the 
standard normalization d?\x) is singular in the Chebyshev case. 

2. Characteristic Polynomials. 

2.1. Jacobi- Tau method. Substituting l|1.6fl into (|1.4|) . the Jacobi- Tau approx- 
imation can be written explicitly in terms of the yet undertermined constants tq, t% 
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and the eigenvalue [i = 1/A, as 

(2.1) u n (x)=T J2» k D™Pi">V(x)+T 1 £ ^D^P^Hx). 

k=0 k=0 

The boundary conditions w„(±l) = then yield the characteristic equations 

If] Pi^l 



(2.2) 



To 



5>^ 2fc p^)(i) + p k D 2k p£;P (l) = o, 



fe=0 



r ^ M ^ 2fe p(«'' 3 )(-l)+r 1 £ f&lpM (-1) =0. 

fc=0 fc=0 



Equation l|BT|) shows that P^'^-l) = (-1)™P^' q) (1) so the 2nd equation above 
can be rewritten at x — 1 by flipping the indices a and (3, 



(2.3) 



r ^ M ^ 2fe p(^)(l)+r 1 £ fti*pM (1 



0, 



fc=0 



r ^ M ^ 2fe p(^)(l)-r 1 ^ /•p 2fc P^J 1 Q) (l) = 

fc=0 fc=0 

This system has a non-trivial solution (tq,ti) ^ (0, 0) if and only if 



0. 



(2.4) £^ 2fe P^)(l) £ ^P^l) 

fe=0 fe=0 

[^] [?] 
+ y k D 2k P ( tf\l)Y j ^ k D 2k P^' a \l)=Q. 

k=0 k=0 

This is the characteristic equation for the eigenvalue fj,. 

2.2. Gegenbauer-Tau method. Gegenbauer polynomials (x) are the class 
of Jacobi polynomials Pn*'^\x) with equal indices a — [3 = 7 — 1/2 (sect. [BT2Tl . The 
Gegenbauer polynomials are even in x for n even and odd for n odd 1, 22.4.2]. Cheby- 
shev and Legendre polynomials are Gegenbauer polynomials with 7 = and 1/2, 
respectively. The symmetry of the differential equation (|1.1|) and of the Gegenbauer 
polynomials allows decoupling of the discrete problem into even and odd solutions. 
This parity reduction leads to simpler residuals and simpler forms for the correspond- 
ing characteristic polynomials. The residual in the parity-separated Gegenbauer case 
contains only one term 



(2.5) 



R n (x)=r XG^(x), 



where Gn^ (x) is the n th Gegenbauer polynomial and n is even for even solutions 
and odd for odd solutions. Substituting (|2.5f) into (|1.4I) provides the Gegenbauer-Tau 
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approximation to in terms of an undertermined constant r and the eigenvalue 
fi = 1/A 



The boundary condition w n (±l) = leads to the characteristic polynomial equation 



since by symmetry Gn'(—1) = (— 1) TI G„ (1) and the two boundary conditions give 
the same equation. 

3. Zeros of characteristic polynomials. 

3.1. Stable polynomials and the Hermite Biehler Theorem. The general 
approach to prove that the eigenvalues are real, negative and distinct is to construct a 
particular stable polynomial p(z) then to use the Hermite Biehler theorem to deduce 
that the polynomials and ^(/-O such that p(z) — ili(z 2 ) + zf^z 2 ) have real, 

negative and distinct zeros that interlace. 

Definition 3.1. A real polynomial, p(z), is a stable polynomial (or a Hurwitz 
polynomial), if all its zeros lie in the open left half-plane, i.e. their real part is strictly 
less than zero, $tz < 0. 

Definition 3.2. Let f2i(/x) and ^(/-t) be two real polynomials of degree n and 
n — 1 (or n) respectively, then f2i(/i) and f^O-*) form a positive pair if: (a) the roots 
fix, . . . , n n of fix and fj,[ ■ ■ ■ , fin—i ( or lAi ' ' ' > tKi) °/^2 are real, negative and distinct; 

(b) the roots strictly interlace (or alternate) as follows: 

fix < n[ < ■ ■ ■ < n' n _x < Mn < (or Hi < fix < - ■ • < /jf n < /z„ < 0); 

(c) the highest coefficients of flx(p) and ^(m) are of like sign. 

We will use the following theorems about positive pairs, p. 198], |17l Sec. 2]: 

Lemma 3.3. Any nontrivial real linear combination of two polynomials that form 
a positive pair has real roots. 

Lemma 3.4. Let P(p-) and Q(fJ-) be real standard polynomials (i.e. the leading 
coefficient is positive) with only non-positive zeros. Then P(/i) interlaces (or alter- 
nates) Q(n) (in the sense of defjnition \3. l A but not strictly) if and only if for all A > 
both Q(fi) + AP(n) and Q(p) + AfiP(n) have only non positive zeros. 

The polynomials P and Q of theorem l3.4l do not necessarily form a positive pair 
since they are allowed to have common and/or multiple roots. We call this set of 
polynomials a quasi-positive pair. 

Lemma 3.5. 7, Lemma 3.4] If flx(p), ^(/-O and 0i(/i), 62(m) are two positive 
pairs then the zeros of H(fi) — f2i(/i)02(/x) + £li(\i)<d\(\i) are real, negative and 
distinct. 

Stability fdefmition l3.1f) is very important in temporal discretizations and matrix 
theory y] as well as in analysis (e.g. |10| and references therein). Stable polynomials 
can surface as characteristic polynomials of a numerical method applied on a differ- 
ential equation. A necessary and sufficient condition for a polynomial to be stable is 
given by the Routh-Hurwitz theorem (see, for example, |121 §40], |131 §23]). Other im- 
portant characterizations of stable polynomials are the Routh-Hurwitz criterion and 



(2.6) 




- 



(2.7) 
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the total positivity of a Hurwitz matrix jlU) . although these will not be used here. 
The characterization of stable polynomials that will be most useful here is given by 
the Hermite-Biehler Theorem |14l p. 197].|1(J). 

Theorem 3.6 (Hermite-Biehler). The polynomial with real coefficients p(z) = 
Q,i{z 2 ) + zQ,2{z 2 ) is stable, if and only i/f2i(/z) and f^A*) form a positive pair. 

The Hermite Biehler theorem states that the even and odd parts of stable poly- 
nomials form positive pairs. This supplies us with a very strong tool to prove reality 
and negativity of the roots of certain polynomials. 

Gottlieb and Lustman 7 used the Hermite Biehler theorem to prove that the 
spectrum of the Chebyshev collocation operator for the heat equation is real, negative 
and distinct for a variety of homogeneous boundary conditions. The basic strategy is 
to show that the characteristic polynomial for that method are the even or odd parts 
of a stable polynomial. Our results extend their strategy to a class of Jacobi and 
Gegenbauer Tau methods that includes Chebyshev and Legendre Tau and Galerkin 
formulations. Although the general approach is similar to that of Gottlieb [S] and 
Gottlieb and Lustman 0, the extension is technically non-trivial and there are dif- 
ferences and some corrections. The key steps in 7 is to prove that the polynomials 
(4. 11), (4. 12)] are stable. To do so, Gottlieb and Lustman derive a first order 
differential equation for those polynomials then transform that ODE into an inho- 
mogeneous one-way wave equation (4.13)] and call on the results [§1 (3. 18), (3. 20)] 
to deduce stability. This is not quite correct since the eigenvalue \i here is complex 
hence WN(x,t) in 7, (4.13)] is also complex while Gottlieb implicitly assumes reality 
oiv N (x,t) and R N (x,t) in (3.18), (3. 20)]. 

The proof for (4.11)] can be fixed and generalized as done in (theorem 13. 71 
below) where we deduce stability of the polynomials (|3.1(l below without going back 
to a one-way wave equation. Here, that proof is further generalized (appendix 0) to 
the polynomials l|3.2|l and (|3.3|l below in order to prove our results about the Jacobi 
Tau method for the 2nd order operator. Our proof follows Gottlieb's ideas to derive 
the results |BJ (3.18), (3.20)] although we do not use Gauss integration. 

The proof for (4.12)] does not appear to be correct however and we have not 
succeeded in obtaining a corrected proof. Gottlieb and Lustman do not provide a proof 
of stability for that polynomial (4.12), they state only that a 'similar argument holds'. 
Gottlieb |S] likewise suggests that the proof of stability for 6] (3.11b)] implies stability 
for [SJ (3.11a)] but this is not evident since T\{t) and T2(t) are distinct functions of 
time that are fully determined by their respective solution procedure. Gottlieb also 
suggests that vx(x,t) 6, (3.8)] is directly related to UN(x,t) |HJ (3.2)] by relation |HJ 
(3.8)]. It is true that Tn(x u ) can be eliminated for n = 0, . . . , N — 1 as suggested in 
the derivation of (3.8)], since (1 - x)DT N {x) = 2(-l) Ar - 1 J2k=o(- l ) kc k lT k{x) as 
can be deduced from |BJ (3. 2), (3. 3)]. However this does not imply that the resulting 
dk coefficients '0, (3.8)] deduced from the afc's that solve [HI (3.6)] are the same dk's 
as those that solve [HI (3.10)]. 

Hence, it appears that there is currently no proof of the stability of [HI (3. 6), (3. 11a)] 
and (4.12)], therefore invalidating Gottlieb and Lustman's proof that the eigenval- 
ues of the Chebyshev collocation operator are real, negative and distinct . 

3.2. Important Stable Polynomials and Positive Pairs. Here we prove 
stability of certain real polynomials whose even and odd parts are directly related to 
the characteristic polynomials derived in section [21 for the Jacobi Tau method. 
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Theorem 3.7. Let Pn"'^\x) denote the Jacobi polynomial of degree n, where 
n>2. If — 1 < a < 1 and [3 > —1, then the zeros of the polynomial 

n / r lk \ 

(3-1) ^(M):=E(^^ ) W M* 

k=0 ^ ' x=1 

lie in the left half-plane; that is, <& n {t>) is a stable polynomial. The proof of this the- 
orem is in 0] together with a discussion of its relation to zeros of Bessel polynomials. 
The next two theorems give two generalizations of the above result that are needed 
for this paper. 

Theorem 3.8. Let Pn a '^\x) denote the Jacobi polynomial of degree n, with 
n > 3, then the polynomial 

71 / \ n— 1 / \ 

(3-2) *M ■= E S + (^^V)) , k 

/==o yax k=o yax ' *=i 

is stable for every A > w/ien — 1 < a < and fJ > — 1. 

Theorem 3.9. Let Pn a '^\x) denote the Jacobi polynomial of degree n, with 
n > 3, i/ien t/ie polynomial 

n / rlk \ n — 1 / jfc \ 

(3-3) := E (^^Wj^AVE (SS^*)) _/ 

«s stable for every A > iw/ien — 1 < a < 1 cmd /3 > — 1. 

The proofs of these theorems are technical and they are given in appendix^ Our 
next theorem combines all the above theorems to get an important result. 

Theorem 3.10. Let n > 3 ; then the polynomials 

(3.4) fl^O*) := Y^^D^P^Xl), OfcfO*) == eV^^CO 

fc=0 fc=0 

/orm a positive pair if — 1 < a, /3 < or < a, /3 < 1. 

Remark 1. It was shown in ^ that these polynomials have real negative and 
distinct roots for — 1 < a < 1 and — 1 < /?. TVie important addition of theorem \S.10l 
is that the roots of these polynomials interlace as was conjectured in Oy. 

Proof. Applying the Hermite-Biehler Theorem to the stable polynomials of the- 
orem 13. 71 for a given n and also for n — 1 > 2, proves that the polynomials n^"'^\p) 
and Q^i\fi) have real negative and distinct roots for — 1 < a, fJ < 1 (Notice that 
we can interchange a and 0). Applying the Hermite-Biehler Theorem to the stable 
polynomials of theorems 13.81 and 13.91 shows that the polynomials 

[f] 

(3.5) J2u k D™PfHl)+A £ M fc i? 2 M-f ) (l) = ^ ) (M)+^fc 1 ) (M) 1 

fc=0 fc=0 

[f] Pi 1 ] 

(3.6) E^ jD2fep ^ ) ( 1 ) + ^ E ^M-f (i) = ^ Q,ffl (/0 + ^£f (m) 

have real negative and distinct roots for all ^4 > and — 1 < a, (3 < 0. These results 
provide sufficient information to apply lemma[!0]to deduce that the set of polynomials 
form quasi-positive pair for — 1 < a, (3 < 0. 
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To show that these polynomials form a positive pair recall that both Cl^"'^\fj,) 
and have real, negative and distinct roots by theorem 13. 71 Thus it remains 

to show that they have no common roots. To do so, assume that fi is a common root. 
Then PiGtz) = flt^ip) + ASl!^ (p) = and P 2 (/i) = fi„ a " 3) '(/z) + jlpQ^-f (/*) = 0. 
Since £T2 n a ' /3 ^(^) 7^ and £)n^p(/x) 7^ (both f2„ and fi n _i do not have a double 
root), set A = — DQ, ? a e) ^ or A = — Dn " (a 1 whichever one is positive (one of 

the two must be since /i is negative). But this will imply that DPi(fi) = and 
DP 2 {p) = respectively, a contradiction since Pi((J.) and P2(/z) have simple zeros. 
For the second range of parameters, replace n with n + 1 in theorems 13.81 and 13.91 and 
apply the Hermite-Biehler Theorem. This gives that the polynomials 



(3.7) Y, f i k D 2k+1 P^ + P(l) + A » k D 2k+1 Pi a 'V(l), 

k=Q k=Q 



(3.8) J2^ k D 2k+1 Pr[+i + £ ii k D 2k+1 P^(l) 

k=0 k=0 



also have real negative and distinct roots. Using l|B.7(l these polynomials transform 
to 

[f] l 2 ^] 

(3.9) g / /^ 2fe J R("+ 1 -' 3 + 1 )(l) + A' £ M fe D 2fc P^+ 1 ^ +1) (l), 

fc=0 fc=0 

[f] I 2 ? 1 ] 

(3.10) Y, f i k D 2k Pi a+1 -' 3+1 \l) + A^ £ /D^P^+^+^Cl). 

The new polynomials have real negative and distinct roots for all A' > 0, where 
A' = rl A, and thus the second range < a. B < 1 follows from the first one 

by a simple change of variables. □ 

3.3. Eigenvalues of the Gegenbauer and Jacobi Tau methods. The pre- 
vious subsection provides all necessary information needed for proving reality and 
negativity of the eigenvalues. First consider the Gegenbauer case. 

Theorem 3.11. The eigenvalues of the Gegenbauer Tau discretization of the sec- 
ond order operator with Dirichlet boundary conditions, vroblem \l.l\ are real negative 
and distinct for —1/2 < 7 < 5/2. Also, characteristic polynomials given by successive 
order (i.e. n — 1 and n ) approximations interlace. 

Proof. The eigenvalues are the roots of (|2.7(l . The theorem follows directly from 
theorem I3.1UI since a = /3 = 7 — 1/2 in the Gegenbauer case and the two ranges for 
the indices (a, (3) merge into the single range —1/2 < 7 < 5/2. □ 

Remark 2. The results of theorem XS. 1 1\ are sharp in the sense that well- conditioned 
numerical calculations (sect. Yh*!fy give some complex conjugate pairs of eigenvalues for 
7 > 5/2 and Gegenbauer integration \B.10\) diverges in general for 7 < —1/2. 

Remark 3. As pointed out in the introduction, the Galerkin method with weight 
function W a ^(x) for vroblem \l.l\ is equivalent to the Tau method with weight function 
W a +i ! /3+i(x). For the Gegenbauer case the Galerkin method with weight function 
W-y(x) is equivalent to the Tau method with weight function Wy+i(x) . Thus, a direct 
consequence of ' theorem Vd '. 1 1\ is that the eigenvalues of the Gegenbauer Galerkin method 
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are real negative and distinct for — 1/2 < 7 < 3/2. Again, characteristic polynomials 
given by successive order approximations interlace. 

Remark 4. Theorem VJ. 1 1\ includes the Chebyshev and Legendre polynomials. For 
7 = 0, we have Gn\x) — Tn ^ , where T n {x) denotes the n th Chebyshev polynomial of 
the first kind \14\ p. 19]. Thus the theorem implies that the Chebyshev-Tau method has 
real negative and distinct eigenvalues with interlacing characteristic polynomials given 
by successive order approximations. For 7 = \ we have G„ (x) = P n {x), where 

P n (x) is the n th Legendre polynomial and for 7 = 1 we have that Gn\x) = 
where U n {x) the n th Chebyshev polynomial of second kind. Therefore, the same result 
holds for both the Legendre Tau and the Chebyshev Tau of 2nd kind method. Fur- 
thermore, in the Galerkin case (see remark\^j the Galerkin Chebyshev, the Galerkin 
Legendre and the Galerkin Chebyshev of the 2nd kind methods have real, negative and 
distinct eigenvalues as well. 
For the Jacobi case we have 

Theorem 3.12. The eigenvalues of the Jacobi Tau discretization of the second 
order operator with Dirichlet boundary conditions, vroblem \l.l\ are real negative and 
distinct if —1 < a, (3 < or < a, (3 < 1. 

Proof. By theorem 13.101 the polynomials (Cln"'^ (fi) , $£—i(p>)) form a positive 
pair for —1 < a, (3 < and < a, (3 < 1. Interchanging the indices a, (3, the same 
result holds for (f2^f '"^(/x), fi^"' Application of theorem 13.51 to these sets of 
positive pairs gives that the polynomial 

(3.11) 4^(m) = n^MJfeV) + <# ,a) (M)n£?V) 

has real negative and distinct roots. Equation 1)2. 4JI shows that BlC'^\^) is the 
characteristic polynomial for the Jacobi Tau method. □ 

This paper focuses on the second order problem with Dirichlet boundary condi- 
tions. Naturally, questions arise about the equivalent results for different boundary 
conditions. The next remarks gives some answer to that question. 

Remark 5. Consider vroblem \l.l\ with Neumann boundary conditions i.e. Xu — 
D 2 u — with Du(±l) = 0. The Jacobi Tau method gives Xu — D 2 u = tqP^ 1 (x) + 
TiPj l °^i\x). Notice that X — is an eigenvalue since u = constant is a solution. For 

X ^ 0, differentiate the last equation to get XDu— D 3 u — ToDPn a '^\x)+TiDP^^\x) 
with Du(±l) — 0. Set now v — Du and with the use of \B. 7| ) the equation transforms 
to Xv - D 2 v = r^P^'^ix) + T[P 7 [ a _ + 2 1J:>+1) '(x) with v(±l) = 0. This is the Jacobi 
Tau approximation of second kind i.e. (a, (3) — ► (a + 1, j3 + 1) for vroblem \l.l\ and by 
theorem VS.l'A it has real negative and distinct eigenvalues for —1 < a, [3 < 0. 

Remark 6. Consider vroblem li.il with boundary conditions u(— 1) = and 
Du(l) = 0. Using eauation \2.1\ and following the ideas of subsection \2.1\ we get that 
the characteristic polynomial in this case is Bn (m) = ^n—i^^ (l^)^^2'^ + \l J ')'^ 
k n ^^l"\^)^n-^i'^ +1 \p), where k n = i(n + a + f3 + 1). Since the polynomials 
(Q,n' a \lJ')i ^nl C i) an d {kn^n-\^ + {^i k n -i^n-2^ + ) form positive pairs for — 1 < 
a, (3 < then by theorem \3.5l the roots of B^^ are real negative and distinct for 
-1< a, [3 < 0. 

An alternative implementation is to rescale the domain to [—1, 0] with the Neuman 
boundary condition Du(0) = and to use only the even Gegenbauer polynomials with 
a Gegenbauer-Tau approach on the entire domain [—1,1] and boundary conditions 
it(±l) = since such polynomials automatically satisfy Du(0) = 0. 
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4. Characteristic polynomial recurrences. The previous section shows that 
the characteristic polynomials given by successive order approximations of the Gegen- 
bauer Tau method have real negative and distinct roots that interlace strictly, pro- 
vided that — 1/2 < 7 < 5/2. Reality of the roots as well as the interlacing property is 
an important characteristic of orthogonal polynomials of successive order ^S],^^ p. 
16]. These properties are direct consequences of the three term recurrence relations 
satisfied by orthogonal polynomials |15|.|14|. Here, we derive recurrence relations for 
the characteristic polynomials of the Gegenbauer Tau method and show that they 
consist of three term recurrences plus a constant term in general. The constant term 
vanishes when 7 = 1/2 or 3/2. 

From H2.7fl , the characteristic polynomials for Gegenbauer- Tau approximations to 
the even and odd modes of (|l.lf> are, respectively, 



(4.1) 



Pm(fJ-) 



k=0 



D 2k G™(l), and q^) := £ 



/ D 2k G 



(7) f1 \ 



k=0 



theorem 13. 1 II states that these polynomials have real, negative and distinct zeros for 
— 1/2 < 7 < 5/2 and that the zeros of p m {y) and q m {lj) interlace, as do the zeros of 
<Zm-i(A t ) an d Pm(lJ-)- Recurrence relations for these characteristic polynomials follow 
directly from recurrences for 2nd derivatives of Gegenbauer polynomials which, using 
(IB. 161) twice, read 



(4.2) 



Lr — 

G« = 



D 2 G% ) 
2(7+l) ! 



G\ 



(7) _ 



D 2 G 



(7) 



4(7+l)( 7 + 2) : 



D 2 G 



(7) 



D 2 G 



(7) 



D 2 G 



(7) 



4(7 + n+ l)(7 + n) 2(7 + n + l)(7 + n- 1) 4(7 + n)(j + n - 



1)' 

4.1. Recurrences for even modes. Substituting l|4.2|) with n = 2m into the 

characteristic polynomials p m (/i) defined in (|4.1J) . with Po(a*) = 1 since Gq 7 ^ := 1, 
yields the recurrence 



(4.3) 



MPo(m) ; 
MPi(i") = 

mp™(m) 

Pin- 



2(7 



M) 



-K, 



(7) 



4( 7 + 2)(7 + 3) 2(7+l)( 7 + 3) 



K. 



(7) 



l(m) 



Pm-l(fj) 



4(7+n + l)(7 + n) 2(7+71+ l)(7+n- 1) 4(7+n)( 7 +n-l) 
where, using (|_B. 171) . 



(4.4) 



and 



K, 



K^> := 



(4.5) K™ 



-(7) ._ 

- — 


G 2 7) (l) 
2(7+1) 


27- 
4(7- 


f 1 

M)' 






•(7) ._ 


Gi 7) (l) 








(27 2 +7- 7)(l + 2 7 ) 


2 - — 


4(7 + 2) (7- 


f3) 


2(7 


+ l)(7 + 3) 


48( 7 +l)( 7 + 2) 




G%{1) 






G, (7) (l) 


, G£? 9 (l) 


4(7 H 


-n+l)(7 + 


n) 


2(7 ■+ 


- n + 1)(7 + n 


- 1) 4(7 + n)(7 + n 
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for n > 3. From 1|B.17|) . the latter expression reduces to 

U 6) KM = G (7) (1) = (27-l)(27-3) (fry + n - 3\ 

{ ' n n(n 2 -l)(n + 2) n - 2[ ' n(n 2 - l)(n 2 - 4) V n-3 /' 

for n > 3, and obeys the recurrence 

1 ' n+2 (n + 4)(n + 3) " 4 720 

4.2. Recurrences for odd modes. Substituting (|4.2|) with n = 2m+l into the 
characteristic polynomials for the odd modes q m {^) defined in (|4.1|) . with <7o(m) = 1; 
gives 



(4.8) 



W0(M) = 4( 7+ l)( 7 + 2)-^ ' 
m m {p) = 

gm+lQf) gm(^) 9m-l(M) _ ^( 7 ) 

. 4( 7 +n+l)(7+n) 2( 7 +n + l)(7+n-l) 4(7+n)( 7 +n-l) 



where, using llRTTl) . #f j := [4( 7 + 1)( 7 + 2)]- 1 G^ 7) (l) = [12( 7 + 2)]- x (2 7 + 1) and 
K^ as in l|4.5|l but here with n = 2m + 1. The recurrence (|4.7|) applies here also but 
starting now with = (2 7 - 1)(27 - 3)/120. 

In general, l|4.3|) and 1)4. 8|) are three term recurrences plus the constants K„ , , 
These constants vanish for all n > 3 when 7 = 1/2, the Tau-Legendre method, and 
when 7 = 3/2, the Tau-Legendre of the 2nd kind or Galcrkin-Legendre method. In 
those cases, the recurrences have only three terms hence the corresponding p m (/x) and 
q m (n) sequences of polynomials are orthogonal polynomials |141 p. 13 and references 
therein]. The recurrence 14.7fl for also indicates why 7 = 5/2 is a critical value 
in theorem IXTT1 For 7 < 5/2, (27 + n - 1)(2 7 + n - 2) < (n + 4) (71 + 3) and ifi 7) 
decreases with increasing n, while for 7 > 5/2, increases with n. If id 7) = 

for n > 3, the characteristic polynomial sequences p m (p) an d 9m (a*) satisfy a three 
term recurrence, respectively, therefore they are orthogonal and have real roots that 
interlace. The constants 7^ pulls down or pushes up the successive polynomials 
in the sequences with respect to that orthogonal case. For 7 > 5/2 that shift leads to 
the bifurcation from real eigenvalues to complex conjugate pairs. 

5. Numerical Implementation. 

5.1. Matrix formulation of the recurrences. The recurrences (|4.3|) and (|4.8|l 

for the characteristic polynomials can be expressed in the matrix form 

(5-1) m[Po(m),Pi(m))---] = [Po(m)>-PiO)>- • •] M 

where the semi-infinite matrix M is tridiagonal plus one row. The matrix M is purely 
tridiagonal if 7 = 1/2 or 3/2. The roots of the m-th order polynomial Pmifj) are the 
eigenvalues of the m-by-m matrix M(0 : m— 1, : m— 1). A matlab code, buildGI2 .m, 
is provided in appendix IU1 which constructs the matrix GI2=A/(0 : m,0 : m — 1) for 
both the even and odd modes by direct implementation of formulas (|4.3|> and (|4.8|l 
with H4.7fl . This approach provides an effective and well-conditioned technique to 
compute the Gegenbauer-Tau eigenvalues as illustrated in figure |5~T1 which shows the 
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odd mode eigenvalues for two values of to = MG, the total number of modes, for 7 = 0, 
0.5, 1, 1.5, corresponding to Chebyshev-Tau, Legendre-Tau, Chebyshev-Galerkin and 
Legendre-Tau, respectively. The Gegenbauer-Tau approximations involve an expan- 
sion of the solution u n (x) into odd polynomials up to degree n = 2 MG+ 1 and therefore 
up to degree 2001 for MG=1000. That calculation shows that slightly more than 60% of 
the spectrum is captured with close to machine precision (here double precision IEEE 
arithmetic), demonstrating the excellent numerical conditioning of the formulation. 
Comparing the MG=100 and MG=1000 calculations shows that there is a slight balloon- 
ing of the round-off error at the higher truncation level. This might be explained by 
assuming randomness of roundoff errors with a standard deviation growing like \/MG- 
It would be interesting to obtain asymptotic estimates for the high frequency modes, 
k > 0.6 MG. The largest eigenvalues for MG=1000 and 7 = 0, 0.5, 1, and 1.5 are, respec- 
tively, 4.86 x 10 12 , 1.63 x 10 12 , 7.61 x 10 11 and 4.07 x 10 u . So the Legendre-Galerkin 
method can be said to be slightly less stiff than the other methods. These values are 
consistent with estimates that the largest eigenvalues are 0(MG 4 ) |18j . 




Fig. 5.1. Relative error |A — A e |/|A e | for the entire odd mode spectrum for MG=100 (left) and 
MG =1000 (right). Exact eigenvalues are A e = — fc 2 7r 2 , k = 1,. . . ,MG. The relative errors for 7 = 0, 
0.5, 1 and 1.5 are shown but essentially indistinguishable at this scale. 



5.2. Don't Differentiate, Integrate. Our approach so far has been theoretical 
and focused on the basic eigenproblem (|l.lf) . For more general two-point boundary 
value problems, e.g. nonlinear problems, it would not be possible to obtain explicit 
forms such as (|1.4|l for the discrete solution, and the residuals would not be as simple 
as <|1.6f) or (|2.5|l . For more general applications it is necessary to select explicit bases 
for the trial and test functions and to perfom the integrals l|1.5|) or l|1.7|) by Gauss 
integration. 

One classical implementation of the Gegenbauer-Tau method is to expand u n (x) 
in terms of Gegenbauer polynomials u n (x) — ^2 7 i =0 0,1 G l f\x) and to use the n—1 

Gegenbauer polynomials G^]\x), k = 0, . . . , n — 2 as the test functions in lieu of 
p n -2(x) in l|1.5|) . Those n—1 integrals, computed by Gauss quadrature in practice, 
and the two boundary conditions yield the n + 1 equations to determine the n + 1 
coefficients a^. For the even modes of the simple eigenproblem i|l.l|) . this formulation 
consists of the even expansion 

m m 

(5.2) u 2m {x) = Y d aiG < 2i\ x ) D 2 u 2m {x) = Y,aiD 2 G%\x) 

2=0 1=0 
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with the weighted residual equations 

(5.3) J (D 2 u 2m - Xu 2m ) G%> W™dx = 0, k = 0,...,m-l, 

where = (1 — a; 2 ) 7-1 / 2 . Equations (|5.3|l yield a matrix problem Aa = XBa for 

a = [ao, ■ ■ ■ , a m ] T , where the m-by-(m+ 1) matrix B is diagonal plus one zero column, 
from orthogonality of the Gegenbauer polynomials (fir?. 10|) . and the m-by-(m + 1) 
matrix A is upper triangular with A(k,l) = for fc > I since, from 14.2fl . D 2 G2i{x) 
can be expressed in terms of all the Gegenbauer polynomials of even degree less than 
21. The boundary condition u n (l) = YlhLo a l^2l W = allows the elimination of 
one of the coefficients, ao or a m say. This elimination can be expressed in the form 
a = Ca where a is the column vector containing the remaining m coefficients and C 
is an (to+ l)-by-m matrix consisting of the m-by-m identity matrix plus one full row. 
This yields the generalized eigenvalue problem ACa = XBCa. The structure of the 
resulting matrices (AC) and (BC) depends on which coefficient is eliminated. If a m 
is eliminated, then (AC) is full and (BC) is diagonal. If ao is eliminated then (AC) 
is upper triangular and (BC) is zero everywhere except on the first row and the first 
lower diagonal. 

Many other implementations are possible. For instance, one can use a poly- 
nomial expansion that satisfies the boundary conditions a priori, U2 m (x) = (1 — 
x2 ) T^iLo bi ip2i(x), where (fi2i(x) is an even polynomial of degree 21. Picking <f2i(x) = 

(x), the equations l|5.3|l lead to a generalized eigenvalue problem Ab — XBb where 
this m-by-m matrix A is upper triangular and B is tridiagonal. All of these formu- 
lations are mathematically equivalent; in exact arithmetic they would provide the 
same eigenvalues as the matrix M in H5.lt . However, the formulations just mentioned 
use 2nd derivatives of Gegenbauer polynomials and these methods are plagued by 
roundoff errors that grow like m 4 , the fourth power of the number of coefficients as 
illustrated in figure I5~2l [31 116 | . 

There is one formulation that is numerically stable and leads exactly to the tridi- 
agonal plus one row matrix of eqn. I|5.1|l . That formulation consists in expanding not 
U2m(x) but its 2nd derivative D 2 U2m(x) in terms of Gegenbauer polynomials: 



m— 1 



(5.4) D 2 u 2m (x)= ^cG^i) 



U2m(x) 



1=0 



1=0 



c a 2 G^(x) 



f3x, 



where I 2 denotes double integration. That double integration is easily expressed in 
terms of Gegenbauer polynomials by double integration of the recurrence formulas 
(|4.2I) which gives 
(5.5) 

' ^ T^h) + °° + ** T2G " " 4( 7+ fh + 2) + * + ** 

J2 G (7) ^ _ ^2 p 

2 4( 7 + 2)( 7 + 3) 2(7+l)( 7 + 3) 

Mi) ^(7) 

j2^(7) _ ^n+2 « n-2 



4(7 + n+ l)(7 + n) 2( 7 + n + l)(7 + n- 1) 4( 7 + n)(-y + n - 1) 
+ a n + f3 n x. 
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The constants of integration a n and f3 n can be defined arbitrarily since the a + fix 
terms have been included in l|5.4|l . so let a n = j3 n = for all n. For the even mode 
expansion considered in this section, we have (3 = in (|5.4|l . so only a survives as 
the lone constant of integration. That constant is determined from the boundary 

condition u n {l) = 0, which for Q reads YTJq 1 c i j2G 2i\ l ) + a = °- From <E3> 
with a n = f3 n = 0, one finds that 

m — 1 

(5.6) a = -^ Cl K^ 

1=0 

with the constants K$ as in l|4.5|l and 14. 4|) . Substituting (|5.4|) with (|5.5() into (|5.3[) 
and using orthogonality of the Gegenbauer polynomials (|B.10|I yields an eigenvalue 
problem Ac = XBc where the m-by-rn matrix B is tridiagonal plus one top row and 
the m-by-m matrix A is diagonal with A(k, k) = J^^G^jJ^W^dx > 0. The system 
can thus be rescaled to the form 

(5.7) c = XMc 

where the matrix M = A~ X B is the tridiagonal plus one top row matrix in (|5.1|1 that 
was obtained from the characteristic polynomial recurrences (recall that /i = 1/A and 
that B is tridiagonal if 7 — 1/2 or 3/2) . That matrix which consists of the coefficients 
in (|4.3() or 1|5.5|) (with a n — fi n = 0) together with the constants —K^ that modify the 
first row and impose the boundary condition can now be interpreted as the chopped 
double Gegenbauer Integration operator with Dirichlct boundary conditions. That is 

if f(x) = IXo 1 fl G 2l\ x ) then 9 = M+ f where M+ = M (° : m,Q : m - 1) and 
/ = [/o, ft, . . . , f m ~i] T provides the m + 1 even Gegenbauer coefficients of the double 
integral of f(x) that vanishes at x = ±1. Note that 14.3fl and l|5.7|l provide direct 
interpretations for the left and right eigenvectors of M, respectively. The problem for 
the odd modes is entirely analogous and does not need to be repeated here since all 
the details are available in l|4.8|l and in the matlab code in appendix El which provides 
GI2= M + . The numerical performance of two differentiation approaches based on 
(|5.2I) . and of the integration approach (|5.4() equivalent to (|5.1() . are shown in figure 
15.21 which displays the relative error for the first even mode eigenvalue as a function 
of m = MG for 7 = 0. The integration formulation (|5.4|) was proposed by Greengard 
p. 1077] precisely for the purpose of controlling roundoff errors. This procedure is 
essentially equivalent to the commonly used reformulation suggested in |SJ p. 120], 
§5.1.2]. 

The Legendre Galerkin (i.e. Gegenbauer Tau with 7 = 3/2) integration im- 
plementation corresponds to Ierley's expansion in associated Legendre polynomials 
|llj . For (|1 . ip . and restricting to even modes, Ierley's expansion consists of u n {x) = 
IXo 1 9i (t-x^Gf/^ix), where G ( n /2) (x) cx P^ 1] {x) cx DP n+1 (x) and P n (x) is the 
Legendre polynomial of degree n (appendix lB|l . Now the derivative of eqn. I|B.2(1 for 
a = [3 = gives D 2 ((1 - x 2 )DP n+1 ) = (n + l)(n + 2)DP n+1 , and, since DP n+1 cx 

G ( n /2 \x), Ierley' s expansion satisfies D 2 u n {x) = 9l ( 2l + l )( 21 + 2 ) G 2i /2 \ x ) 

and corresponds to an expansion of the 2nd derivative of u n (x) in terms of Gegenbauer 
polynomials of index 7 = 3/2, a special case of the integration approach l|5.4|) . Ier- 
ley's test functions (I — x 2 )G^J 2 \x) vanish at x — ±1, so his equations are l|1.7J) with 
a = (3 = corresponding indeed to a Legendre Galerkin approach (or Gegenbauer 
Tau with 7 = 3/2). This yields an eigenvalue problem of the form Ag = XBg where 
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Fig. 5.2. Relative error A — A e |/|A e for the first even eigenvalue as a function of MG = m for 
the Chebyshev-Tau method, 7 = 0. The exact eigenvalue A e = — 7r 2 /4. The dashed line indicates 
MG 4 scaling of roundoff errors. Three implementations are shown, the differentiation approach 
15. H6 with a m eliminated (top curve) and with oq eliminated (middle curve), and the integration 
approach |5.^| ). The latter is well- conditioned with errors staying at the level of machine precision 
10~ 15 . The gaps in that curve occur where the approximate eigenvalue is indistinguishable from the 
numerical value for 7r 2 /4. 

A is diagonal and B is tridiagonal, where the coefficients g have been renormalized 
so that B is also symmetric. 

6. Conclusions. It has been shown that the eigenvalues of the Jacobi Tau 
method for the second derivative operator with Dirichlct boundary conditions are 
real, negative and distinct for ranges of the Jacobi indices a and 0. These ranges 
include Tau methods with Chebyshev and Legendre polynomials of the 1st and 2nd 
kinds. Chebyshev and Legendre Galerkin formulations are included as well but col- 
location methods are not. Although our work owes much to earlier work by Gottlieb 
and Lustman d , we have raised doubts about the validity of their proof for the 
Chebyshev collocation operator. 

Special emphasis has been placed on the symmetric case of the Gegenbauer Tau 
method where the range of parameters included in the theorems can be extended and 
characteristic polynomials given by successive order approximations interlace. The 
interlacing is between g m _i(/i) and p m {p) in (14. 1ft - and between p m (/u) and q m (/J>), 
not between p m (n) andp m +i(/i) or between q m (n) and g m +i(//), although we believe 
the latter hold as well [3] conjecture 3]. Proving such interlacings could allow a proof 
for the spectrum of the Gegenbauer collocation operator since the parity-reduced 
residual in that case reads xDGn^ (x) which can be written as a linear combination 
of GL 7+1) (x) and C^\x), from itTTTl and (H/Hf . 

The characteristic polynomials for Gegenbauer Tau approximations have been 
shown to satisfy three term recurrences plus a constant term that vanishes for the 
case of the Legendre Tau and Galerkin methods. Hence for those two particular cases 
the characteristic polynomials are orthogonal, and their roots interlace. A well con- 
ditioned matlab code that computes the roots of the characteristic polynomials for 
general Gegenbauer parameter 7 is provided in appendix El In section 15.21 several 
mathematically equivalent numerical formulations are discussed. The theoretical and 
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practical superiority of the integration method, which is numerically stable, is em- 
phasized. In a forthcoming paper we apply similar methods to the simplified Stokes 
eigenvalue problem D 4 u = XD 2 u with u(±l) = Du(±l) = and rigorously identify 
classes of spectral methods that are free of spurious eigenvalues. 

Acknowledgments. The authors thank Jue Wang for several helpful calcula- 
tions in the early stages of this work. 

Appendix A. Proof of Theorem |3"T81 and Theorem 13.91 

Proof. (Theorem 13 .8|) For A = the theorem reduces to theorem 13.71 Fix now 
A > but otherwise arbitrary. Let 



(A.l) f n (x; »)=J2 ^ k D k P^ {x)+AY J n h D*pW\x) 

with fj, such that /„(l;/i) = and D := d/dx. Then /„(x;/x) satisfies the following 
differential equation 

(A.2) ( /n _pKfl_AP£f>)= M ^. 

Multiplying by ^"j 1 0-+ x )t integrating from —1 to 1 in the Jacobi norm and adding 
the conjugate we obtain: 

(A.3) J 1 + x)W a ,pdx J 1 + ^) (1 + x)P^{x)W a ,pdx 

(l+x)W a: fjdx. 



L 


dfn 




dx 



For the first term of equation l|A.3[) . integration by parts yields 
(A.4) J ^^(l + x )W a ,pdx = - J ^ IM 2 ^fy {fi + 1 - a - ((3 + 1 + a)x) dx. 

For /3 > — 1 and a < the factor [3 + 1 — a — ((3 + 1 + a)x is nonnegative for all 
x e [—1,1]. The second term of ljA.31) can be expanded as 

(A ' 5) j\ + W) (1 + x ^ )W ^ dx = 2 j\ xDPt^Pt^W^dx = 

2B n . 1 t xP^P^W a ,,dx = 2i? "- lfl1 - 1 f (p^Yw^dx = 2B -^- l f 
7-i 03,n-i 7-i v ' 



O3 n-l 



where we have used expression ()B.8(1 , the Jacobi recurrence relation (|B.5|) and orthog- 
onality of Pn"'^\x) to all polynomials of degree less than n with respect to the Jacobi 
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weight W a ^(x). Similarly, the third term of equation i|A.3|) can be calculated as 
(A-6) f ( d -k + d -i) (l + x)P^W a , pdx 



_i \ dx dx 

2 I (1 + x)DP^P^ ] W a ,pdx + fa + /!*) J xD 2 P^P { ^ ] W a ,pdx 
+ 2A C xDP r [ a JpP,[ a JpW a , dx 



/■ 2 /* 2 

[Pi-P) W a , dx + 2B n ^j x[P n ^ W a ,pdx 

+ (/i + /i*)P„_ 1 S n _ 2 y xP£$P<*pW a ,f,dx + 2AB n _ 2 J xP^P^W^dx 



= 2B n _ 1 h^ 1 -2^^B n _ 1 h^ 1 + {^*)^ 

<23,n-l 03 jn _2 a3 jn _2 



a,/3 



Substituting these expressions back into equation (|A.3|I yields 



(A.7) 



(1-*) 



(/3 + l- a-(/3 + l + a)x)da; + 



-l<2l,n— 1 , Q 



0^3, n-1 



+ 2A ^1- 

= (m + m*) 



»3,n-l 



Pn-1 + A 



03,n-2 



",/3 

n-2 I 







[/'. 


da; 



(1 + x)W a . p dx + A^^P n _ 1 P n _ 2 / l «f 1 

0*3, n-2 



Since 1 > °^' T1 - the left-hand side is positive and < 0, which ensures stability. 
□ 

Proof. fTheorem l3.9|) For A = the theorem reduces again to theorem 13.71 Fix 
now A > 0. Let 



JV 



iV-l 



(A.8) f n (x;n) = ^n k D k P^>(x)+A^ E ^P^-fix). 

with / n (l;/i) = 0. Then f n (x;fi) satisfies the differential equation 
(A.9) i^-^'-A^f 



/i 



dx 



Multiplying by f*(x,fj,) \i_Jy integrating from —1 to 1 and adding the conjugate 
yields 

1 f 1 ^ 1 + x ) P (^) {x)Wa ,edx 



(A.10) (- + ^)[ Ifnfp^rW^dx-l I f 
M / J -i 0--x) A* J-i 



(1-x) 



~ f fn^\P n a - 0) (x)W a ^dx-A^ I J 

M J-i 



1 ~< 1 + *>p£f(*)W^dx 



-i (1-^ 

1 d|/»| 2 (l + ^) 
. x (1 — x) 



Wa^dx. 
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Integration by parts on the first term gives 



(A,l) / 



— (l+x)W aJ ,dx=-Jjf n \ 



(/3 - a + 2 - (/? + dx. 



For /3 > — 1 and a < 1 the factor /3 — a + 2 — (/3 + a)x is nonnegative for all a; € [—1,1]. 
For the other terms on the left hand side of (|A.10(1 . recall that /«(!; fx) = so write 



(A.12) 

then 

(A.13) 

Also 
(A.14) 

rl 

-'71—1 



1 £ ± 4fnPt' 0) (x)W a , dx= f c* n _ 1 xPt?{x)P^\x)W a , dx 
1 I 1 — XJ J-l 



Ct3,n-1 J-i v 7 a3 jn _i 



1 i^/^f^w^dx 

1 I 1 - X J 



+ x)P^P(x)P^'(x)W a , dx + I c* n _ 2 xP^(x)P^(x)W a ,pdx 



>(«,/8) 



1 (i^-f^V^dr-cj;.! 



Ct3,n— 1 J-l 
1 



-n-2" 



i£f (*)) Wa^dx 
P^f(x))V a)/3 dx 



1 c n-l 


^ 02,n-l 




<J3,n-l. 



_* Ql,n-2 \ , a,p 



-n-2" 



03,n-2 



Explicit values of c„_i and c„_2 follow from equation (|A.12|I 



(A.15) f n = {l-x)Y J c k Pi a ' P \o 



k=Q 



>-iP£?\x) - c^xP^f (x) - c„_2xP<«f (x) + 0(n - 2) 



-Cn-V 



Pi a ^(x)+(c n . 1 



Q2,n-1 
Ct3,n-1 



C„_2" 



Ol,n-2 
; 

Ct3 : n-2 



Pi a J\x)+0(n-2). 



0-S,n-l 

Now from equation (|A.8J) 

/„ =P^(x) + nDP^\x) + An 2 P^\x) + 0(n - 2) 
=P^\x) + B n _ xl iP^\x) + Afx 2 P^J\x) + 0(n - 2). 



(A.16) 



Comparing these two expressions for /„ gives 
(A.17) 



C n -1 = 



a 3,n-l 



Cn-2 



a3,n-2 
0>l,n-2 



02,n-l 03,n-l „ ,2 

B n -lH - Afl 

1l,n-l dl,n-l > 
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Substituting all these results back into (|A.10|> yields 

,(!+*) 



(A.18) 



1 1 

— I 



f n \ 2} h-^W a . dx + h^ 



(1-x) 



n-1 



2|,,|2,a,/3 _ 



or after rearranging some of the terms 



2 (a-0 + 2-(a + 0)x) 
(1-x) 2 W ^ dx 



(A.19) 

H 



' (1-x) 
1 , 2 (0-a + 2-(0 + a)x) 

\Jn\ 



2A 



W a ^dx-2AB n _M 2 K-i 



(1-x) 2 

The right hand side is negative so this implies that 5i(/i) < 0. □ 
Appendix B. Jacobi and Gegenbauer polynomials. 

B.l. Jacobi Polynomials. The Jacobi polynomials Pn (x) are suitably stan- 
dardized orthogonal polynomials on the interval (—1, 1), with weight function W a ,[3 — 
(1 — x) a (1 + x)P . The class of Jacobi polynomials Pi°'^ (x) includes Gegenbauer (Ul- 
traspherical) polynomials when a = (3, Chebyshev polynomials when a = = — 1/2 
and Legendre polynomials when a = = 0. 

Definition B.l. The Jacobi polynomial, Pi a '^\x), of degree n, can be defined 

by 



(B.D e»> M :4Etr)(l!f) 



(x-l) n - k (x + l) K , a, > -1, 



where the binomial coefficient (?) = (a)(a — 1) •• • (a — k + l)/k\. Jacobi polynomi- 
als are the most general class of polynomial solutions of a singular Sturm-Liouville 
problem on the interval — 1 < x < 1 and this is directly related to their excellent 
approximation properties §9.2.2, §9.6.1]. The Jacobi polynomial P^*'^\x) satisfies 
the differential equation 



dx 



(B.2) -f- ( (1 - x) a+1 {l + xf +1 ^y ) = n(n + a + + 1)(1 - x) a (l + xfy. 



i d 



dx 



Jacobi polynomials IjB.ljl are orthogonal with respect to the weight W a ,p{x) = (1 
x) a (l + xf 



(B.3) 

where 

(B.4) 



{l-xT{l + x)?P^P^dx=\ °^ 



h?>f> 



2 a+/3+i T(n + a + l)T(n + /3+1) 
2n + a + 0+1 n\T(n + a + + 1) 
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Orthogonal polynomials satisfy a three term recurrence relation, for the Jacobi poly- 
nomials this reads 

(B.5) 2(n+l)(n + a + (3 + l)(2n + a + (3)P^\x) = 

((2n + a + + \){a 2 - /3 2 ) + (2n + a + /3) 3 x) P^\x) 

-2(n + o)(n + /3)(2n + a + /3 + 2)p£f >(*). 

where (2n + a + (3) 3 = (2n + a + /3)(2n + a + /3 + l)(2n + a + /3 + 2). To ease the 
notation in calculations we write the recurrence relation in the form 

(B.6) a hn P^\x) = (a 2 ,„ + a 3 , n x)Pt^(x) - <n, n P<*?\x). 

Two other useful relations involving derivatives of Jacobi polynomials [5] are 

(B.7) T X P ^ ){X) = \^ + * + P+ l)P^' 0+l \x). 

and 

(B-8) ^ P «+i /3)(x) = B ^ P) ( X ) +Pn-x{x) 

with B n = ( 2 "+ a +^+iK^~^^+^+ 1 ) anc j p n _ 1 (x) a polynomial of degree n — 1. 

B.2. Gegenbauer Polynomials. The Gegenbauer (a.k.a. Ultraspherical) poly- 
nomials ci 7 "*^), 7 > —1/2, of degree n are the Jacobi polynomials with a = (3 = 
7 — 1/2, up to normalization 22.5.20]. They are symmetric (even for n even and 
odd for n odd) orthogonal polynomials with weight function W(x) = (1 — cc 2 ) 7 ". 
Since the standard normalization 22.3.4], is singular for the Chebyshev case 7 = 0, 
we use a non-standard normalization that includes the Chebyshev case but preserves 
the simplicity of the Gegenbauer recurrences. Set 

(B.9) :=1, GW(x):=^, n > 1. 

27 

We refer to these non-standard Gegenbauer polynomials as ns-Gegenbauer for short. 
The ns-Gegenbauer polynomials satisfy the orthogonality relationship 

(B.10) [ l (l-x 2 y- 1 / 2 G^G^dx= I °' m * n > 

J -1 1 K, m = n, 

where P 22.2.3], 

mm 7r2- 1 - 27 r(n + 2 7 ) 
(B.llj K - —, ■ n . 2( . ■ 

The derivative recurrence formula 1|B.7|) for ns-Gegenbauer polynomials reads 
(B.12) ^ G (-r) i=2(7+1)G (7 + D j 

(for Cn^ this is formula [21 A.57]), and the three-term recurrence takes the simple 
form 

(B.13) (n + l)G^l 1 ^2(n + 1 )xG ( ^ ) -(n-l + 2 1 )G ( ^l 1 , n > 2, , 
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with 
(B.14) 



G[, 7) (x) = 1, G™(x) = x, G£> = (7 + l)s 



-.(7) 



(7) 



Differentiating the recurrence (|B.13|I with respect to x and subtracting from the cor- 
responding recurrence for 7 + 1 using l|B.12|) . yields ^| 22.7.23] 



(B.15) 



(„ + 7 )GW = (7 + l) 



G (7H 



G 



(7+1) 



n > 3. 



Combined with IB. 121) . this leads to the important derivative recurrence between 
ns-Gegenbauer polynomials of same index 7 

G$\x) = ^'(4 2(l+ 7 )G^(x) = 



_d_ r 

dx 



G 



(7) 

71+1 



G 



(7) 



(B.16) 2(n + 7)GW 

Evaluating the Gegenbauer polynomial at x = 1 we find ^ 22.4.2], 

(B.17) G(7)(l) = J-CW(l) = J-( 27 + n - 1 

27 27 \ n 

where f^" 1 ) = (2 7 + n - 1)(2 7 + n - 2) • • • (2 7 )/n! = ^g^. 

Gegenbauer polynomials correspond to Chebyshev polynomials of the 1st kind, 
T n (x), when 7 = 0, to Legendre P n (x) for 7=1/2 and to Chebyshev of the 2nd kind, 
U n (x), for 7=1. For the non standard normalization, 

(B.18) G (°) {X) = IM 7 G ^\x)=P n (x), GW(,) = ^, 

n l 

Appendix C. Matlab code for Gegenbauer- Tau Double Integration. 

function GI2=buildGI2(MG,g, ip) 

7, buildGI2 produces the Gegenbauer-Tau double integration operator GI2 with 

7« Dirichlet boundary conditions u(+/-l)=0 for even (ip=0) or odd (ip=l) solutions. 

7. 

7. GI2 = buildGI2(MG,g,ip) yields the (MG+l)-by-MG tridiagonal + 1 row matrix GI2. 
7o (2*MG+ip) is the degree of the polynomial expansion, g is the Gegenbauer index 
7o g=0 is Chebyshev-Tau, g=l/2 is Legendre-Tau, g=l is Chebyshev Galerkin, 
7« g=3/2 is Legendre-Galerkin. g must be greater than -1/2. 
7. 

°/o EXAMPLE: Cheb-Galerkin odd mode eigenvalues compared to exact values: 

7, MG=20; GI2=buildGI2(MG,l,l) ; M=GI2(1 : end-1 , : ) ; eCG=sort (1 . /abs (eig(M) ) ) ; 

7. k=[l:MG]; semilogy(k,k. ~2*pi~2,k,eCG, 'o' ) 

7. 

7. Fabian Waleffe & Marios Charalambides , 2005, 2006 
n=2*(l:MG-l)+ip; 

dm=l./(4*(g+n+l) . * (g+n) ) ;d0=-l . / (2* (g+n+1) . * (g+n-1) ) ; dp=l . / (4* (g+n) . *(g+n-l) ) ; 



T=diag(dm(l:MG-2) ,-l)+diag(d0)+diag(dp(2:MG-l) ,1) ; 7. Tridiagonal part 
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7« K_n by recurrence (minus sign included) 

if (MG>2), Kn=zeros(l,MG-2) ; K3=(2*g-l)*(3-2*g)/120; 

if (ip==0), Kn(l)=(4*g~2-l)*(3-2*g)/720; 7„m=2, n=4, Kn(m)=K_{2m+2} 

elseif (ip==l), Kn(l)=K3*(2*g+2)*(2*g+l)/42; 7.m=2, n=5, Kn(m)=K_{2m+3} 
else error (' ip must be or 1'), end 
for m=2:MG-2; 

n=2*m+ip ; Kn (m) =Kn (m-1 ) * (2*g+n-l) * (2*g+n-2) / ( (n+4) * (n+3) ) ; 
end, end 

7 1st row and 1st column 

if (ip==0), M00=-(2*g+l)/(4*g+4); 

M01= (7-g-2*g~2) * (l+2*g) / (48* (2+g) * (1+g) ) ; M10=l/ (2*g+2) ; 
elseif (ip==l), M00=-(2*g+l)/(12*g+24) ; 

M01=l/ (4* (g+3) * (g+2) ) + K3 ; M10=l/ (4* (g+1) * (g+2) ) ; 

end 

rl=[M00, M01, Kn] ; cl=[M10; zeros(MG-2, 1)] ; re=[zeros(l,MG-l) ,dm(end)] ; 
GI2=[rl; cl,T; re] ; 
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